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Abstract 

Worm algorithms have been very successful with the simulation of sigma 
models with fixed length spins which result from scalar field theories in 
the limit of infinite quartic coupling A. Here we investigate closer their 
algorithmic efficiency at finite and even vanishing A for the one component 
model in dimensions D = 2, 3, 4. 
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1 Introduction 



In reference [1] it was shown how spin systems of the Ising and XY model type 
can be Monte Carlo simulated with greatly reduced or even eliminated critical 
slowing down by what the authors call worm algorithms. They can be seen as 
stochastically sampling strong coupling or hopping parameter expansion graphs 
of these models to in principle arbitrary order instead of the original spin or field 
configurations. The method has been generalized in several directions, including 
sigma models of the O(A^) type, and a summary with further references can be 
found in [2]. All these systems can be regarded as lineai0 sigma or 0^ models, 
where in the limit of infinite self coupling A the length fluctuations of the spins are 
frozen to — 1 dimensional spheres. In this paper it is our goal to systematically 
investigate for = 1 the dynamical behavior of worm algorithms away from the 
Ising limit at finite A and even in the Gaussian limit A = in Euclidean space 
dimensions D = 2,3,4. Results for the Ising limit A = oo may be inspected for 
comparison in |3] and The paper is organized as follows. In Section 2 we 
discuss the reformulation of 0^ theory and of a few observables by the all-order 
strong coupling contributions. This is followed by Section 3 detailing our update 
algorithm and presenting our results on critical slowing down. In Section 4 we 
summarize our conclusions and in Appendix |X] we compile data for A = 1/2 that 
can be useful as a reference for future studies. Our report here is a summary of a 
diploma thesis available under [5] with a lot more details. 



2 (Re-) formulation of the model 
2.1 Partition function 

We use lattice units (a = 1) and start from the standard lattice formulation of 
scalar field theory with the partition function 



j^i=(.y) mm 



where the measure at each site z is given by 



d^^x{ci>)f{<^>) 



oo 
oo 



f{<p) 



(2) 



and thus includes the coupling A ^ 0. Our hypercubic lattice is wrapped on a 
D dimensional torus with L sites in each direction, and the sum is over links / 



^ Spins live in the linear space but, of course, are still nonlinearly coupled in general. 
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corresponding to nearest neighbor pairs of sites {xy). Below we shall need the 
moments 

CA(n) = j rf/iA(0)0". (3) 

In particular for the A values simulated in this study the nonvanishing even mo- 
ments are 



Co(2n) = r(n + l/2)/r(l/2) = 2-"(2n-l)!!, 
ci/2(2n) = 2"/2r(n/2 + l/4)/r(l/4), 

Coo(2n) = 1, 

while odd moments vanish due to Z(2) symmetry. 

The fundamental two point correlation function is regarded as a ratio 



(4) 
(5) 
(6) 



G{u-v) = {ct>{u)ct){v)) 



Z{u, v) 



(7) 



with the 'partition function with insertions' as numerator 



Z(u, v) 



YldfixiHz)) 



By expanding in /5 the Boltzmann factor on each link and then integrating over 
independently on each site using the same quantity is given by 



Z{u,v) = 



n 



x,v ) 



(9) 



Here the sum is over a link field that is independently summed over k{l) 
0, 1, . . . , oo and the divergence 



dk{x) = J2 K^) 

l,dlBx 



(10) 



counts the sum of k{l) over all links surrounding a site x. We have a well defined 
non-negative weight for any k configuration. Nonzero contributions arise if dk is 
even everywhere except the sites u, v if they do not coincide, where dk must be 
odd. Finally we introduce the ensembl^ 



k,u,v 



n 



J]^ cx{dk{x) + 5x,u + 5x,v) 



(11) 



weight depending on the locations u,v could be included as in [J, but we have not yet 
explored this generalization. 
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which will be simulated by the worm algorithm. In this way expectation values 



k'a,v 



n 



X.l 



0[k-u,v] (12) 



become accessible to Monte Carlo estimation. 
2.2 Observables 

The fundamental two point function is now given by 



G{x) = (0(x)0(O)) 
where we have introduced 

rx{2n) = 



cxi2n) 



cx{2n + 2) 



(13) 



(14) 



which is only required and defined for even arguments. For universal quantities 
only ratios of two point functions at differing separations are of interest and then 
the denominator is not required. It is however related to the susceptibility 



(15) 



The second moment (renormalized) mass m can be defined and then measured by 



cos(27ra;^/L)G'(x) 



m 



with 



m2 + 4sin^(7r/L) 
fm{x) = ^^cos{2nx^/L) 



{{fm{u-v)))^m (16) 



(17) 



Note that in our symmetric setup each of the D directions contributes equally and 
we average for better statistics. 

An estimator for the energy (without the measure part) is given by 

p^^Ai 7 1 {{5^,.rxmu))Y.ikm ..o^ 

A second estimator for E is given by the nearest neighbor correlation using f fT3|) . 
As in the Ising limit |1] the error is much larger in this case. 



4 



3 Algorithm and dynamical results 



In the first Subsection we describe details of tlie worm algoritlim implemented by 
us. Tlie code was written in C and lias run on standard dual quad-core PCs. The 
data reported below correspond to a few core-months total runtime dominated by 
the largest lattices of sizes 128'^ and 32^. 



3.1 Specification of the update scheme 

The update algorithm used by us is a generalization of the one described in [4J. 
We define two types of micro-steps that each obey detailed balance with respect 
to the ensemble (fTTl) : 

• /: With equal probability we pick one of the 2D neighbors u' of the present 
configuration's u and denote by / the link in between. Again with equal 
probability we propose one of the two moves k{l) — )■ k{l) ± 1 accompanied 
in either case by the change u — ?■ u' . The proposal is accepted with the 
respective Metropolis probabilities min(l, qi±) where we take 

^i+ = ^^^^'\\9k{u') + 5^.,,) (19) 

and 

qi_ = ^r^{dk{u) + 5u,.-l). (20) 

• ir. We only act iiu = v holds, and even then only with probability po = 1/2. 
Then a randomly chosen new location u' = v' (at unchanged k) is proposed 
and accepted with probability min(l,g//), 

_ rx{dk{u)) 
r\[ak[v!)) 

Swapping the roles of u and f , step / can also be applied to move v and we call 
these possibilities now /„ and /„. As an iteration we define the L^/2 fold repetition 
of the sequence J„ II J„ II which is similar to a sweep of a standard local algorithm 
as far as CPU work is concerned. 



3.2 Monte Carlo dynamics in the Gaussian Umit 

For the Gaussian case A = all our observables can be computed exactly by 
straight forward Fourier expansion. We therefore do not list any mean values in 
this case but we have monitored that our results were correct within errors. The 
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dynamics of the 'worm' algorithm remains of interest at A = and is presumably 
representative for other small bare couplings. To investigate critical slowing down 
we keep constant the extension in physical units mL while we increase L. In the 
free case this leads to choosing 



D + — . 
2 



(22) 
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Figure 1: Integrated autocorrelation times for the energy E (left panel) and the 
mass m (right panel) for the Gaussian model in dimensions D = 2,3,4 with 
mL = 4. The dashed lines are fits of the form Tint oc L^, values are in the text. 



All autocorrelation times were determined as discussed in [6] and are given 
in units of 'iterations' defined above. Measurements are taken and pre-averaged 
during each iteration and then stored for off-line analysis. In Fig. [T] we see log- 
log plots of integrated autocorrelation times for the observables E and m. Note 
that the estimator f|T8|) refers to a ratio of primary Monte Carlo estimated mean 
values and we refer to [6J for the definition of Tint for such derived quantities. The 
dashed lines are fits with dynamical exponents z. From top to bottom we have 
determined z = 1.62(2), 0.97(3), 0.51(2) (left plot) and z = 1.86(7), 0.99(2), 0.29(2) 
(right plot). The quoted errors are purely statistical. We consider these fits which 
have have acceptable over the range shown as mere parameterizations of our 
data in the range where most simulations work. We have not embarked on the 
difficult assessment of systematic errors with regard to truly asymptotic dynamical 
behavior. In summary we see here distinct critical slowing down, weaker than 
for standard local methods but inferior to worm (and also cluster) simulations 
in the Ising limit. There is a pronounced tendency of decorrelation improving 
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with increasing dimension. The susceptibihty x was also investigated and behaves 
similarly to m. The ultraviolet quantity G(0) measured via (fT3|) on the other 
hand typically exhibits shorter autocorrelations. We also have a large set of data 
in smaller volumes mL = 1. Autocorrelations are larger for this more critical series 
but the overall qualitative behavior is quite similar. Further details can be found 
in 0. 



3.3 Monte Carlo dynamics at intermediate coupling 

We now turn to simulations at A = 1/2 which are otherwise organized similarly to 
those of the previous subsection. As the values of our observables are non-trivial 
we collect them in Appendix |A] for this case. The dynamical results are given in 
Fig. El 
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Figure 2: Integrated autocorrelation times for energy E (left panel) and mass m 
(right panel) for A = 1/2 in dimensions D = 2,3,4 with mL = 4. The dashed 
line is a fit for Tint linear in InL, while the dotted lines just connect data points to 
guide the eye. 



We have chosen semilogarithmic plots as the most natural appearing common 
format here. Two term fits linear in InL do however not quite achieve good 
values except for E in D = 4 (dashed line). What is clearly visible is that the 
worm algorithm works much better in this interacting case then in the Gaussian 
limit. Although precise fits require more than the first power of InL (or larger 
L to assume the asymptotic form) the general behavior looks close to an only 
logarithmic growth of autocorrelations. The tendency with D is mild and the 
ordering differs for the two observables displayed. Also at this coupling we have 
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looked at further observables and have also simulated mL = 1 with qualitatively 
similar outcomes. 

4 Conclusions 

We have found that a worm algorithm is applicable to 0^ theory all the way to 
the Gaussian limit in dimensions between two and four. Although still exhibit- 
ing less critical slowing down than standard local methods, the free limit A = 
is the most difficult case for the worm algorithm. For the intermediate coupling 
A = 1/2 the situation is already similar to the infinite coupling limit with regard 
to decorrelation. While at zero coupling the efficiency shows pronounced improve- 
ment with growing dimension, this dependence is much weaker at A = 1/2. In [7] 
a very dramatic improvement was realized in the Ising limit by the construction 
of a special estimator for the connected four point function without the need to 
perform numerical cancellations. We have so-far not succeeded in generalizing this 
construction to finite A in a similarly efficient way. Some such efforts are however 
reported in [5]. 

Acknowledgments. We thank Martin Hasenbusch for discussions. Financial 
support of the DFG via SFB transregio 9 is acknowledged. 

A Results of simulations at A = 1/2 

In this appendix we report on some mean values of observables defined in (fT6|) . 
(1181) . fllSp together with /3 values that have been determined by tuning to mL = 4 
within errors. The data are summarized in tables [H 121 El The last column holds 
the numbers of iterations performed. We note that the mass m has a practically 
constant and mostly sub per mille error at - or rescaled to - a constant iteration 
number. This confirms the very mild or even absent slowing down. 



L 


/3 


mL 


E 


X 


its/10^ 


8 


0.576950 


4.0036(20) 


0.26184(13) 


6.566(5) 


24 


12 


0.615670 


3.9982(22) 


0.32214(13) 


13.381(11) 


24 


16 


0.634350 


4.0008(23) 


0.36241(12) 


22.274(20) 


24 


22 


0.649270 


4.0028(25) 


0.40381(11) 


39.334(38) 


24 


32 


0.661390 


3.9991(26) 


0.44622(9) 


76.855(80) 


24 


64 


0.674330 


4.0067(29) 


0.50603(7) 


262.87(30) 


24 


128 


0.680679 


3.9936(35) 


0.54516(5) 


898.3(1.2) 


20 



Table 1: Resuhs at D = 2. 
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L 


/3 


rriL 


E 


X 


its/10^ 


8 


0.370240 


3.9994(26) 


0.14561(7) 


10.610(12) 


8 


12 


0.383030 


3.9984(27) 


0.16235(5) 


22.915(28) 


8 


16 


0.388310 


3.9957(28) 


0.17178(4) 


39.963(50) 


8 


22 


0.391920 


4.0007(29) 


0.17989(3) 


74.06(10) 


8 


32 


0.394390 


3.9999(30) 


0.18692(2) 


154.08(21) 


8 


64 


0.396400 


3.9998(32) 


0.19457(1) 


600.02(90) 


8 


128 


0.397067 


4.0059(49) 


0.19816(1) 


2336.9(5.3) 


4 



Table 2: Results at D = 3. 



L 


/3 


rriL 


E 


X 


its/lO^ 


4 


0.245490 


4.0021(20) 


0.079119(52) 


4.0522(34) 


8 


8 


0.271670 


4.0012(17) 


0.092157(16) 


14.624(14) 


8 


12 


0.277630 


3.9977(15) 


0.097502(8) 


32.183(23) 


8 


16 


0.279870 


4.0003(15) 


0.099984(5) 


56.669(39) 


8 


24 


0.281560 


3.9972(14) 


0.102197(3) 


126.869(83) 


8 


32 


0.282173 


4.0007(19) 


0.103132(4) 


224.41(20) 


4 



Table 3: Results at D = 4. 
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